{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian Data Analysis, 3rd ed\n",
    "##  Chapter 6, demo 3\n",
    "\n",
    "Posterior predictive checking  \n",
    "Light speed example with a poorly chosen test statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os, sys\n",
    "# add utilities directory to path\n",
    "util_path = os.path.abspath(os.path.join(os.path.pardir, 'utilities_and_data'))\n",
    "if util_path not in sys.path and os.path.exists(util_path):\n",
    "    sys.path.insert(0, util_path)\n",
    "\n",
    "# import from utilities\n",
    "import plot_tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# edit default plot settings\n",
    "plt.rc('font', size=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# data\n",
    "data_path = os.path.abspath(\n",
    "    os.path.join(\n",
    "        os.path.pardir,\n",
    "        'utilities_and_data',\n",
    "        'light.txt'\n",
    "    )\n",
    ")\n",
    "y = np.loadtxt(data_path)\n",
    "# sufficient statistics\n",
    "n = len(y)\n",
    "s2 = np.var(y, ddof=1)  # Here ddof=1 is used to get the sample estimate.\n",
    "my = np.mean(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A second example of replications\n",
    "nsamp = 1000\n",
    "pps = np.random.standard_t(n-1, size=(n,nsamp))*np.sqrt(s2*(1+1/n)) + my\n",
    "# Use the sample variance as a test statistic\n",
    "# This is a poor choice since it corresponds directly to\n",
    "# the variance parameter in the model which has been fitted\n",
    "# to the data.\n",
    "pp = np.var(pps, axis=0, ddof=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAEUCAYAAADeE7fRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FEX+P/D3Jwc5J5CQC0IACS6HHMlyqOjKoYscfhUW\nBEUFRFZ2PQF1BUVAFEFXcRVvRCAiCqI/WIOAosCqiAgE8ECQyJ2EK+S+k/r9UT2hM8wkk5Cr4f16\nnnmS6a6urq7pnv5MVXW3KKVAREREZFUe9V0AIiIiogvBYIaIiIgsjcEMERERWVqlwYyIjBURJSJt\nK0izWEQOVacAInJIRJa6WY5x1VlHQycifYw67lPfZakvIjJTROp1AJexLy42ve9jlMvDIV1r4/Ma\nX+eFrAV1vf+5851iRca+0q8W8x8iIpNrML8LKq+ITBSRv7nIt0rHck3mRZemmmqZeQbA0BrKy5Wx\nAC7KYIYajKHQ+7JdHwAzwBZMcs8MALUWzAAYAqDGghlceHknAjgvAAHwLoCr6zEvugR51UQmSqmk\nmsiHqD4ppRLruwx1SUQ8AUh9l4MuLkqpYwCONbS86OJWI784nXUziUgbEflcRHJF5KSIvCQi9xrN\ny62d5HGbiOwVkRwR2S4i15rmbQLQG8A1xvLKmOaqPIEiMl9EjohIgbH+DSLS3pRGichsEXlSRI6J\nSJ6I/E9EYp3k9zcR2WpsS7qIfCwiLZ2ku1dEdotIvoicFpGFIhLikCZMRJaJSKaRVzyAJq5r97x1\n9BaRr0Qky6ir9SLSyTS/k7Et/3FYbrZRF3823vuKyMsi8rOIZItIqoh8Zq4jI529S6CXiKww1ntC\nRKYa8weISKJRlh9FpJvD8ptE5FsRucVYV4GI/CYiI9zYVi8RmWqkLxCRZGM/8q1kuc9EZIPpvYjI\nKSMPf9P0D0TkR9P7sm4mEZkJ/csVAIrs+53DqjxFZJaIpBif5Wci0sKN7XK7Toz6/d74TDNEZJWI\ntHNIIyIySUT2iUihUZ7XRCTIIZ19n58iIgcBFALo7GSd843P2Nthus34/OdWsn0BIjJXRJKMbUsV\nkU9EJMIhaajxGWQan+2rjp+tiDQTkXjjeCoQkT0icqdDmkgRWWLkUWBsf4KIhJvS+IvI8yJy0Kij\ng6KPfQ9TGnt3281G/Z02XktFpMJj1LRvPCnnvqNmmuZXeNwaaW4UkS3G55xtfJ7TjXmLAYwBEGXK\n/1AF5fESkWeMz8D+ffStGN+rFZVXRHqIyEo59724T0SeExE/U/6HALQCcIdp+cXGvPO6hkTkYdHf\n73kiclb0d/zQaublJSKPi8ivxradEpF14vDdRZcYpVSFL+juHQWgbQVpFgM4ZHrfCEASdEQ9BsAg\nAJ8COGzk1dqU9pAx/UcAwwHcBCARQDqAJkaajgB2AtgN4Crj1bGC8iwAcALAPQCug+4+eBHAVaY0\nCsBRAN9BN9+OBLAPwBkAIaZ0/zDSvmdsx0gAewEcBGAzpZsLoAjASwD6A7gbwHEAPwDwNKX7BkAm\ngAcA3Gjke9RYR59KPovBAIoBrAZwi/HaAuAsgGhTuvsAlAIYaLzvB6AEwGRTmsbQTbi3QQeKQwF8\naeQV6eTz/x3AUwBuAPC2Me15AD8ZedwE4FdjWxqZlt8EINX4jO82tiHBKF9fU7qZAJTD9n4EIAfA\ndGO9Dxr7xSeV1NNkALkAfIz3XY315QPob0qXDOB5h31xsfF/C6N+FIBrYOx3xrzWxvRDAJYBGAi9\nn58GsMmNY8rdOhlgfG5fArgZwCgABwCcAhBlSvecUZ7XjH1qEoBs6H3Nw2GfP25MH2bkHwHdnVa2\n/0EfbwrACIdyTzDKeFkF29YIep/MMfaXv0If1wsAtHeyT80yPtunjG192pRXAID9xvbea9TzB8ay\n95rSfWmkuwP6eL8VwFswvmegW6C/gT62JwK4HsCTxv7wkikfez0cBDAf+jh+EEAegCWVfKZXGcsu\nwrnvqBbuHrcA2gAoMLZvAPQxOwHG/gkgBsAaACdN+cdVUJ4njX3gYejj+/8APA3gZjfKOwzANOhj\nujf090kqgI9M+ccBSAGwzrR8jLNj2fhciqGP477Q36NTANxT1byMaSuN/F406moIgHkwHTt8XXqv\nyhNUL5i511imp2maQAcjzoKZswCCTdO6G+lGmaZtAvCtWxsF/AxgXiVpFPTJJ8A0rTV0QPKM8T4Q\nQAaA9xyWvQz6V+1E03IlAKY7pLvGWM8Q4/1fjfe3OaRbC/eCmQMAvnKYFmRsx38cpq+GDug6QZ/A\n1gGQCvL2BOAPIAvAJCef/3TTNC/oL9UimE5s0CdcBaC3w+emUD6Q9ATwG4BvTNPKfWkB+Iux3GiH\nct5hTI+tYFvizOWAPoHtgT7pzTGmtTfSDHDYFxc7lgmAl0P+rY3pmxymP2pMb17J5+hunWyHPuF7\nmaZdZtT7PON9CPRJcLHDOu401nGzwz6fDMDPIW0fx/3PKKPjvrYTwLpKtm2c43qdpLHvU087TE8A\nsN/0/gHHchnTNxj7n6fxPhvAQxWs7y4jn+scpj8JfRyHO9TDEod0r0EHPi6PH1P9Plud4xY64FMA\ngirIfzGAYxWVwaEuP61OeR3SCPTxfid0INvU4XhZ6mSZmSh/LL8GYGcl63E3r35GuV1+3nxdmq/a\nGth4FYAjSqlt9glKKQXgExfpv1dKnTW9/8n4e15Xjpt+BDBWRJ4Qke6ixwY487lSKsdUxkMAtuLc\ngLOrob90PjCaNr1ExAu69eE36F+BgA5SPJyk+wE6OLjOlF8Jzq+HjyrbIBG5HPrXmeM6cgF8b1qH\n3Tjok9526C+jMcZnYM5zhIj8ICLp0L90cqADuHY431r7P0qpYugv6P1KqYOmNL8Zf6Mdlj2qlNpq\nWr4EwMcAeorDlUImA6BPNCsdtvcLY77j9prtBpCGc4Mb+wH42niZpxUB+LaCfCrzucP7quy3FdaJ\niAQA+DOA5UZ929MdhG5N7G1Mugq6NcTxisCPoD/T3g7T1yml8two3xsA+hr7HUSkB3SQ+HYly/UH\nkKqU+q8b61jj8P4nlK+76wAcV0ptcki3FEAYdAsSoI/3x4yujM4i4jgOaAB0K9gWJ/uSN3QdVlYu\nH+hWrCqpwnG7C3p//EhEhpu7yKrpRwCDRHcrXisijapQ5iCjSy4JOlAuAvA+dGBzeTXLEiu6+/IG\nMXX1VkN/6GBmwQXkQReh2gpmmkH/cnJ0wkX6NPMbpVSB8W+FYyMq8CD0l+446APppOjxIY4HkbPy\nnAAQZfxv/0LZAH1Am1+dATR1SHfASTqbKV0zAGeVUkVulMORfR0LnazjJtM6AABKqTPQX8o+AD5U\nSpVbh4j8H4Dl0F1mowBcCaAHdJO+s3o/6/C+0MU0OFneVT03gj4pORNuzM9B+W2171dNXSwHpVQp\ngM3QJ2NP6BPGRuPVTfRYkr4AflRKZbvKxw1pDu+rst9WVifB0CePFCfpUqFbZGD6Wy6dEQCdMc2H\ns3QV+H/GeiYY7/8B3arzWSXLNYVuCXSHs/rzMb0Pgevtt88HdNfvfwH8C7oF7riITDcFyuHQYzIc\njxv7jy3HfelCPldHbh23SqkD0F2EHtCBQ6rocXqOwai7noMe83UzjC42EVkkIqFuLLsI+vN+FfqH\nWg8A9xvzqlMH8QD+Cf0dsx5Amoh8Kk7GTrqhKYA0NwNyuoTUyNVMTqTg3K8msyr/sqkO4wQ1FcBU\nEWkF3YQ7F/pk+3gl5YnAuS/jM8bfsQB+cZI2yyFdf5x/gjfPTwEQLCLeDgGNO/Viz2MqdHDlqND8\nRkRuAPB36JaZ+0RkqVJquynJbQAOKKXGmpbxxvknv5rgqp4LoYMnZ85AN+3/xcX85ErWuRG6T/1a\n6NamzdDdEbnQrRV9UHkrQ22qrE78oH+BRjpJF4lzJ9w007SyfdT49d8U55+YFdyglCoSkXeh950X\noPeXl8ytRC6chu7arAlpcN5KGGmaD6XUSeiT7f2iB0ePgR4fcgrAm9D70kEArgadH6qh8jrj9nGr\nlNoIYKOI+EB3Uc8CsEZEWiulTldlpcb3y/MAnheRSOjAaR50V/JIV8uJHoB9C4CZSqlXTNPPGyhe\nhbIo6GPtbREJhv6efAn6x9SVVczuNIAQEfFjQENmtdUysxVASxHpaZ9gNP0Ou4A8C6C/4KtEKXVY\nKfUSdFOx45fsIKM5317G1tBNzt8bk7ZAByxtlVLbnbz2Gem+hO5Pbukinb0r5nvosRGO9XCbG5uy\nD/pL9woX69hj2o5Q6F9DnwPoBT2gepmIBJry84fuhjC7yyhfTYsWkbKmfKO15FYA24xWFGfWQf8K\nbOxieysLZr6GbuV4CkCiUirdOBH/D3pQZCh0wFMR+y/yKu93bqiwTozuzx0AbjV3kxrBeS/oMS2A\nPtYKcf4+NBL6x8omVN/b0FfafQzdYuJO0/4XACKNlr8LtRlACxG5xmH6KOgWul8dF1BK7VNKPQH9\no8J+vK+D7vrMdrEvVSlQqEAhzt9X3D5uTdtQoJT6GsAL0IOgLzNmVfc7MFUp9S50MGX+DnRWXh/o\n7wDH1uOxTrKucnmUUmeVUssBrHAoi7t5fQHdYnlR3LCSak5VWmYGiEiqw7QMpdSXTtIuhm4B+VRE\nnoT+hTQeuukc0Cf+qvoV+lfiSOgrpbJMwUQ5IvI9dLPzT9C/xntDX9GyxCFpHoAvROTf0Afx09BX\nGr0MAEqpTBF5DMDrIhIGPW4kA7obqjf0ANBlSqkkEXkewGvGL8PN0K0K0dDNtO8qpTYqpb4UkW+h\nf6GEQg/uHAk3fskqpZSI3A9gtdH/vQL6V0oE9MntiFJqnpH8PegD/m7jF/Yo6IBmPvTVM4D+gh8i\nIi9DDxbsjnNXC9W0EwCWi8gM6H3hnwD+ZPx1Sim1SUQ+hB4zMw+6S6AUevDtIACPK6X2V7D8LyJy\nEvrKlX+bZtlbbAqgx55UxH6yfERE1gIocWjduhDu1MlT0F2FCSLyBnQL09PQ++BLAKCUShORl6Bb\nIXOgA9gOAJ6FHg/kOP7DbUqp4yLyX+gr3T5TSh11Y7Gl0C2CH4rIHOhxYzboLpT/KKV+q2hhB4uh\nA0/798gx6AHgfwUwQSlVIiKNoU/SH0CP2SqCblkIxrnxVR9A7/dfGXW1GzrQjYHuhhmilMqtQrlc\n+RXAYBFZBx1MJSulkt05bkXkH9DdoZ9Dj8kLhW7NSYa+oMGef4iI/BO6xTVfKfUTnBCR1cZ27jTK\nEgc9dsjcGumqvFuh9/kUo6zjcK7r3XF7/yIiN0F3/Z02xh06luUd6B+F30MHoX+C/uH0hSmZW3kp\npTaKyCcA5olINPSPFm+j7tY4GV9FlwrHEcGOL5y78sDZ62cjzWKYrmYypsVAH5h50F/Wr0AHOAr6\n17Y93SE4H8WuoJs67e8jjfyy4ORKEodln4c+eWdAj7n4CQ6j3408ZgN4AvpLMh+6b/m8q2SgT54b\noQOdXOgg5D04XB4OfYBuNdaZDT0e5TUYlzwaacIAfGhsRzp0C8otcONqJmP5q6GDj7NGmQ9BD/a8\n2pj/APRJ/68Oy9mvbhlpvPeAPuElG9u0GfoL7xDKX9Fj//zbOuS3CQ5Xl+HcVT7jHdNBnzR+hg4i\n9tnLYUo3E+dfgukBfTLbbWxrhvH/C+Z9qIK6Wo7zr1iyX+l03v7jZNs9AbwO/QVcirIW8/O305je\nx53P0d06MdIOgD4J5BnbvxpAO4c0An059j7oX9spRrmDHNK5utrGZbkB3G7MG1xZfZuWCYQOIA+b\nyrMS564acrVPOdsHmkGPITlt1NMeAHea5vtAn6B/gT7mMqHHyY1yyMfXyP83I580I91MGFeLmerh\nBodl7eVtXcl2XwPdmpaP87+/KjturzY+26NG+VKgW8TamfIIgP7uOGvkf6iCsjwC/V10xth39hnb\n6l1ZeaH377XQ31Enob/DBjvuI9BXBH4D/f2hcO62BuU+R+huv01GXgXQXX4vw7R/upuXMc0L+kq0\n/TjXLfs5HI4Lvi6tlxg7R50QkQQAHZRSMXW2UtdlUQBmK6Wm1XdZLmaib27opZS6trK0lwor1YmI\nfAB90mujXHcJEhHVq9oaAAzRD0TLhm7FsEGPBxiMCroWiKhhMMbzxEJ3g05mIENEDVmtBTPQzYmT\noO8b4QndzDleKbWwFtdJRDXje+gfI0ug7zlDRNRg1Wk3ExEREVFNq61Ls4mIiIjqBIMZIiIisjQG\nM0RERGRpDGYsQETGiogyvbJEZLeIPGDctt7dfF41Lo+HiLR3yNPVa42z5Rs6EWkuIktE5LRRX8tF\npIlp/kQR+UlcP+jyQtff4Or3QurEKOvYaqyzVuvZYV3RIrJSRDJEJFP083+q/LBaEVlnbO+zDtOH\ni8gnInJYRPJEZJ+IzBERW81tBRFVB4MZa7kV+uZaw6DviDsfwHR3FhSRGOiHx800JqUaedlfDxrT\nn3SYfr+L5RssEbkMun6CoO8Y+0/oO9C+Zkr2NvQNDMfUUjEaVP3WY53Udj0DAEQ/RPZr6JuvjYG+\ngeXl0M86CqhoWYd8boe+W7gzj0I/9f4J6JsZvgldj1/WRbBGRBWo77v28VX5C67vmLoR+pESrpbz\nMf0/H/op0a7SPmiso6OL+RUu31Be0HfD3Qp9t1UxTZ8FfbsAX9O0FwD8UgPr9HUjTb3Vb03UiVH2\nsdVcf43UcyXreBg60GhrmnYZ9PPHJruZRzB0EGq/6/GzDvPDnCwz2kjbrza3jy+++Kr4xV8T1vYj\ngCARCReRmUbTeCcRWS8i2dDPgYHop/DeCWBZBXnFQt/W/LznXblaXkS+MJ7j4pi+s4gUicgd1d6y\n6hsK/STeyUop830HjkA/j6e5adpHADqKSK+qrEBEPESkp4jMMLbfnecf1Wf91nidVLFs1arnKroZ\nwFal1AH7BKUf8Pod9ONC3PE89CNaPnQ2Uynl7AnvPxp/nT27iIjqCIMZa7sM+tdotmnaaujnLN0M\n44GZ0E8CbwL97BNXukJ/kZc4medq+e8AxBknYwBlT0d/A8AWpdQH7m9KjRkHfcO3P0TEy/6Cfl4Q\nUP5J4bugnz8zoLJMRSRUREaJyPvQv95/MNa1C/oZX5Wpz/qtjTqpStkqzFM0LzdeFT3R/QqceyCj\n2S8AOlayLRCRa6FbWe6vLK2D3sbfvVVcjohqUG3eAZhqnqdxErIBGAHgb9BPM87V5xEAwKtKqVcc\nlrsKuil8j7NMjTyvgH66sDOulv8O+pd9HHQ3BqBPCFcZ0+qU6KcS9wXgD/30ZEdF0A/WBAAopUpF\nZDd0eZ3lFwbgPgADAfSADhy/hX6I4udKqV/cLFe91W9N10l1yuZGnr2hu0wrsxn6YZDOhEA/gNFR\nGnT3kUtGHb0N4EWl1HktZxUsFwXdVbdB1dzT1ImoGhjMWMtvpv9LoU+OEx3S/D8nyzUHkKmUKnSR\nbzvopwrvcjHf1fJboU/wVwHYalwZ8wKA15RSzn4l17aO0Cft+6EHu5p9BOCsUqrYYfopAH9ykV9X\nnBuQuxfAFABrlVLOgoKK1Gf91nSdVLdsFeW5AzpYrEyWG2mq418A/OBeCxsAQEQCoVtBiwHcXUvl\nIiI3MZixlqEAjkF/qR9WSuU7SZPiZJov9EBPV2KNv7tdzHe6vFIq2+EX92zoIGuGs0xExMvJibMm\ntTb+fquUKmvlEJEI6C45Z2Mh8qBPZM5sBtAfumVmIPTJK1NENgBYC2CdUuqYG+Wqk/p1obXxt6bq\npLplqyjPbLgO9MqttoJ5Z+G8BcZViw0AwLh0+0kA4wH4mLvNjPdNAGSZuwdFxA/AZwDaAOjt5j5A\nRLWIY2as5Wel1Hal1D4XgQzg/Av/DPSYDFe6Gsu5OtlWtPx3AK4SkT9DX1r8mFIq0z7TGJQ8Q0R+\nBPCKMa27iHwlIttFZJeIjHZIP1NEdorIARGZUEG5HdmDc8dxKXcZ27fYyTIhAE47y0wpVaSU+lIp\nNVkp1QH65PUEAB9jW46KyB4ReaiSctVa/bqhRuvkAspWUZ69obu7Knt9VUFZfoHuynPUEcCvFSzX\nBjqYXAod9NhfgL4U+yyAzvbEIuINYCWA7gAGKaV+qiBvIqojbJm5NPwGoJGItHDxKzIWwKEKTkQV\nLf8t9GXH8QC+U0otdZaBUqoHABi/dBcCGKyUOma83y4iP5jGKzRSSv3ZGJOwU0S+c7Nb5ZDx9wro\nkxtEJBLA4wDeUUolOVnGfv+VShlXx7wO4HUR8YU+CQ8E0KWSRWu9fitwyPhbG3VSlbJVlGdNdDP9\nF8CLItJGKfUHAIhIawDXQHcPurILekyRo43QAc5CAAeM/Dygu3b7AbhJKXXe1VxEVD8YzFwa/mf8\n7QndTeWoK/Sv7Oosb1+uPYA/u1h+oen/XtBdHwmmQcte0L+g7cHMOwCglDouImuhTx4/A2UnqIMA\nnlZKzXRYzw7osS1zRCQfugXlGeiT0aOOhTICqT8BeNFZoUXfiK2iO8geBvAWgNwK0gC1XL91WSdV\nLZs7eSqlsgBc6ADaBQAeALBaRKZBtzo9A+Ao9OBee1l6Q7fwjFNKxSul0gFsclJmQHflmue9Dn3j\nytkAckTEPKD5GLubiOoPg5lLgFLqkIhsA/B/AD41zzN+pYfDdRdIhctDj3coBPCmeUyGkzRlq4S+\ngVpV7jli7jqz38011Uk5lYgMhT6xrYDu1lgK4BmllLOAY7BRdmeDpgEdeH3pRvlcXmVTR/Vbl3VS\n1bJVNc9qUUrliEg/6NsRvA+9n30FYKJSynH/80T1utgHGn+fNF5mT8MCd8cmumjV9137+KqbF/Rd\nhDMA+Nfk8gBegh503NjFcgpAE9P7YCP9ANO0jgACTemfNf5vBn2C7mRKey/0lTHV2g6Hsq0F8H59\nfzYXUr91VSdwcgdgd8rW0OqZL774ujhfHAB86VgKfT+R+y50eRHxF5GrReRf0LeRv08pleFOJkqp\ns9C/1B8T/bDMXwG8Cn3PErsiEdkJfRO5Gar8eJneAF5WzlsV3CYisdDdV09fSD416ELqt87qpKpl\na4D1TEQXIVGqoqsd6WJi9PH/WSn1xoUsDz2uYzWA4wDmKKVer8EyKgDBSo9lqDUiMsBYj9Nb19eH\nuqjfStbvsk6Mz+Vu6JvQuV22hljPRHTxYTBDDUpdBTNUNfZgRim1uL7LQkTkiMEMERERWRrHzBAR\nEZGlMZghIiIiS2MwQ0RERJbGYIaIiIgsjXcAPocjoS91ffrov5s21WcpiGqDVJ6EyLrYMkNERESW\nxmCGiIiILI3BDBEREVkagxkiIiKyNAYzREREZGkMZoiIiMjSGMwQERGRpTGYISIiIktjMENERESW\nxjsA0yXvniXbAQCPpWYBAP5tvHfHwjHda6VMRETkPrbMEBERkaUxmCEiIiJLYzBDRERElsYxM9Tg\n3FOFMStmHL9CRHRpYssMERERWRqDGSIiIrI0BjNERERkaQxmiIiIyNIYzBAREZGl8WomImoQMjMz\ncfLkSRQVFdV3US46hYWFh+q7DEQXoATAt8XFxX/v1q1bobMEDGaIqN5lZmbixIkTiIqKgp+fH0Sk\nvot0sTld3wUgqq7S0lI5fPjwtenp6f8E8IqzNAxm6KJR3fvTUP07efIkoqKi4O/vX99FIaIGxsPD\nQzVv3jw7KytrLFwEMxwzQ0T1rqioCH5+fvVdDCJqoBo1alSklGrsaj6DGSJqENi1RESuGN8PLmMW\nBjNERHXsyJEjCAwMRElJSX0XxS3Z2dnSr1+/tjabLXbgwIFt3FmmZ8+e7ebNmxda22VzZtSoUS0f\ne+yxZjWdtiL79u1rJCLd3B3APmzYsNYPPfRQ8wtdL2kcM0NEVIEBAwagZ8+emDVrVrnpq1evxoQJ\nE3Ds2DF4eVXtq7Rly5bIzs6uyWLWqiVLlgSfOnXKOy0tbZe3t/d58ydPntw8KSnJZ/Xq1QfroXjn\nWbZs2ZHaSFtfevbs2e622247M3nyZA7kdoEtM0REFRgzZgyWLl0KpVS56e+//z7uuOOOKgcyxcXF\nNVm8OnH48GGfNm3a5DsLZBoaK9YvXTgGM0T14J4l26v1oro3ZMgQnDlzBt98803ZtLNnzyIhIQGj\nR48GAKxZswZxcXEICgpCdHQ0Zs6cWZb20KFDEBEsXLgQLVu2RL9+/cqm2U+8ixYtQocOHWCz2dCm\nTRu8/fbbZctv2rQJLVq0wEsvvYTw8HA0a9YMixYtKpufl5eHRx55BK1atULjxo1x7bXXIi8vDwCw\ndetW9OrVCzabLbZdu3YdExISbK62c+fOnb49e/ZsZ7PZYtu2bXvFBx980BgAJk2a1Pzll19utmbN\nmmB/f/+4l19+uVzX0cqVK4Pmz58faZ/frl27jvZ5hw8fbvTnP/+5fUBAQNw111xzeUpKSlnk99VX\nXwXExcW1v5CyAbq75o477mjZu3fvtn5+fnEJCQk2xy6cadOmRYSFhXUJDw/vMm/evFAR6fbzzz/7\n2Je3p01ISLBFRER0mTFjRkRISEjXsLCwLq+88kpTez4fffRR4w4dOnQMDAyMi4yM7DJ58mS3u4m+\n++47v44dO3YICAiIGzx4cJuCgoKy8++pU6c8+/bt2zY4OLhrUFBQbN++fdsmJSV5A8CDDz4YtWPH\njsCpU6e29Pf3jxs9enRLALj77rujIyMjuwQGBsZdccUVHdatWxfoblkuRgxmiIgq4OfnhxEjRiA+\nPr5s2ooVK9C+fXt07doVABAQEID4+Hikp6djzZo1ePPNN7Fq1apy+WzevBl79+7F+vXrz1tHeHg4\nEhISkJmZiUWLFmHSpEnYuXNn2fzU1FRkZGTg+PHjWLhwIe6//36cPXsWAPDoo49ix44d2LJlC9LS\n0vDCCy/Aw8MDx48fx+DBgzFt2jSkp6fvmjt37rE777wzJjk5+bympIKCAhkyZEjbfv36ZZw6dWr3\nvHnzjtx7771tdu/e7fPyyy8nP/jgg6mDBw8+m5ubmzhp0qRyXR3Dhw/PNM/ft2/fr/Z5n376acji\nxYsPnjjRwHwmAAAgAElEQVRxYldRUZHHM888EwEABw8e9B42bNjlU6ZMSbmQstnT/Pe//w158skn\nU7KzsxP79+9frv9u5cqVQW+99Vbk2rVr9yclJf28efNml0ETAJw5c8Y7IyPDMyUlZc9rr712eMqU\nKS1PnTrlCQCBgYGlS5YsOZiRkZG4evXq35csWRL2/vvvN6koPwDIz8+XW2+9te3IkSPPpKWl7Ro+\nfPjZdevWlS1XUlKCMWPGnD5y5MhPhw8f3uPr61s6YcKElgAwf/784926dcueM2fOkdzc3MT4+Pgj\nANCjR4+cXbt2/XL27NnE4cOHp915550xubm5l+woegYzRNTwTJwI9OlTu6+JE90uzpgxY7By5Urk\n5+cDAOLj4zFmzJiy+X369EHnzp3h4eGBLl264Pbbb8fmzZvL5TFz5kwEBAQ4vQR98ODBiImJgYig\nd+/e6N+/f7mWIG9vb0yfPh3e3t4YNGgQAgMDsW/fPpSWluK9997DK6+8gqioKHh6eqJXr17w8fHB\n0qVLMWjQIAwaNAienp4YOnRoZqdOnXI++eST8y5v3bhxY0Bubq7n7NmzU319fdXNN9+c1a9fv/Ql\nS5Y0dUxbFbfffvuZLl26FAQGBqq//e1vaT/99JM/ALz77rtN+/TpkzFy5MiMmijbDTfckN6/f/8c\nT09P+Pv7l+sPXL58ecjIkSNPd+/ePd9ms5XOnj07uaIye3l5qX//+9/JPj4+auTIkRl+fn6le/bs\n8QWAm266Katnz555np6euPLKK/NuueWWtE2bNlUYHNm3obi4WJ566qmTPj4+6u677z7buXPnXPv8\nyMjIkrFjx6bbbLbS4ODg0qeeeipl27ZtFeZ73333pUVGRpZ4e3vj6aefPlFYWCi7d+/2rawsFysO\nACa6AOz6uTRce+21CA0NxapVq9CjRw9s27YNn376adn8H374AVOmTMHPP/+MwsJCFBQU4NZbby2X\nR3R0tMv8165di6effhr79+9HaWkpcnNz0blz57L5TZs2LTc2x9/fH9nZ2Th9+jTy8/MRExNzXp6H\nDx/Gxx9/jM8++wwlJSWxAFBcXCzXXXddlmPao0ePekdGRhZ6enqay1uYnJx8QYNkIiMjyy7t8ff3\nL83NzfUwytZo7dq1wTabrSx4uZCytWjRwuUlRKmpqd7dunXLsb+PiYlxejt8u8aNGxebxwb5+fmV\nZmVleQDA119/HTB16tSo/fv3+xUXF0thYaHHwIEDz1aUn30bwsPDizw8zrUftGjRosD+f1ZWlseE\nCROiN23aFJSZmekFADk5OR7FxcUux2RNnz49YunSpaGnTp3yNtJ7njx58pI9p1+yG05EDdh//lPf\nJTjP6NGjER8fj3379uHGG29ERERE2bxRo0bhgQcewNq1a+Hr64uJEyfi9OnyF564uo9OQUEBhg0b\nhvj4eNxyyy3w9vbGkCFDzhtw7ExoaCh8fX2RlJRU1uVlFx0djbvuugsLFiwAgF0V5RMdHV2Umpra\nqKSkBPag4ejRo40uv/zygoqWM21b5YUtv77CoUOHnvnoo48Ou5G20rJVtP6IiIiiY8eONbK/T0pK\nauQqbWXGjBlz2fjx409u3Ljxd39/fzVu3LjoM2fOVHoejYqKKjp58qR3aWkp7AHN8ePHfS677LIC\nAJg1a1bEgQMHfLdu3bq3ZcuWxVu2bPG75pprOtr3AcftW7duXeBrr70WuW7duv3dunXL8/T0RFBQ\nUKw7+8zFit1MRERuGD16NDZs2IAFCxaU62ICgKysLISEhMDX1xfbtm3DsmXL3M7X3pITFhYGLy8v\nrF27Fl988YVby3p4eGDcuHGYPHkykpOTUVJSgu+//x4FBQW488478dlnn2H9+vUoLi5Gbm6uJCQk\n2OwDS8369OmT4+vrW/rUU09FFhQUSEJCgu3rr79uctddd6W5U46IiIjiY8eONXL3vjn33HPPmQ0b\nNjT55JNPgmq7bCNGjEhbvnx50507d/pmZWV5TJ8+vdr3lMnJyfEMCQkp8ff3Vxs3bvRftWpViDvL\nXX/99Tmenp5q9uzZ4QUFBbJkyZIme/bsKXt2R1ZWlqevr29paGhoyYkTJzxnzJhRbmBxWFhY8R9/\n/FE2RigjI8PTy8tLRUZGFhUVFcmjjz7aLCcnxxOXMAYzRERuaN26NXr16oWcnBzcfPPN5ea98cYb\nmD59Omw2G2bNmoURI0a4na/NZsOrr76KESNGIDg4GMuWLTsv/4q8+OKL6Ny5M3r06IGQkBA8/vjj\nKC0tRXR0NFavXo3nnnsOTZs2jY2Kiury4osvRpSWlp7XROTr66tWrVr1+5dfftk4NDS068MPP9zy\nzTffPBgXF5fvThlGjx6dBgDBwcGxHTt27FBZ+rZt2xatWLHiwNy5c5vVdtlGjBiROX78+JP9+/dv\nFxMT0+nKK6/MMfItdWd5s5deeunInDlzmgcEBMTNmjWr+U033VRpF5N9G5YvX560bNmy0JCQkNgV\nK1aE3Hjjjen2+VOmTDmRn5/vERoaGnvllVd26N+/f4Z5+YkTJ55ISEgIDgoKih07dmz0sGHDMnr3\n7p3ZsWPHztHR0Z19fX1LIyMjK+w+u9jJpdws5YAV0UDU1ziUx+ZMAAD8e+rblaSsPwvHdK/vItSK\nvXv3okOHSs+BVH076rsADcXOnTt9e/TocUV+fv4OK9w3h87ZvXt3aNeuXVs7m8eWGSIiuqjFx8c3\nycvLk1OnTnk++uijLfr27ZvOQObiwmCGiIguagsWLAgLDw/v2rZt286enp5q4cKFDf4RBlQ1vJqJ\niIguat98883v9V0Gql1smSEiIiJLYzBDRERElsZghoiIiCyNwQwRERFZGoMZIiIisjQGM0REdezI\nkSMIDAyEu7f/b8ief/75sKZNm3b19/ePS01NrfSW+q+++mrTbt26tauLslXV5MmTm99yyy2XAcDv\nv//eyN/fP664uLhey1TV+oqKiuq8atWqSp/kfbFhMENEVIEBAwZg+vTp501fvXo1IiMjUZ2TXcuW\nLZGdnQ3zk6CtqKCgQGbMmBGdkJCwPzc3NzEyMrJcdLZv375GItKtqMjlQ60brMsvv7wwNzc30dVT\nq93Vs2fPdvPmzQutoWLVKBHp9vPPP/tUnrLm1NY+wfvMEFGDVNuPtXD30RBjxozBk08+iaeffrrc\nk6/ff/993HHHHajqya64uLjKyzRUx44d8yooKJDu3bvn1XdZAKCoqAi8s++liS0zREQVGDJkCM6c\nOYNvvvmmbNrZs2eRkJCA0aNHAwDWrFmDuLg4BAUFITo6GjNnzixLe+jQIYgIFi5ciJYtW6Jfv35l\n0+ytOosWLUKHDh1gs9nQpk0bvP32ueeDbdq0CS1atMBLL72E8PBwNGvWDIsWLSqbn5eXh0ceeQSt\nWrVC48aNce211yIvT8cWW7duRa9evWCz2WLbtWvXMSEhoaz74dVXX23aokWLzgEBAXFRUVGd33zz\nTadPgM7Ly5Nx48ZFh4eHdwkPD+8ybty46Ly8PNmzZ49Pp06dOgFA48aN46666qo/OS7bp0+fdvb5\n/v7+cRs2bAiwz7v33ntbBAUFxUZFRXVesWJFkH36mTNnPEeMGNEqLCysS3h4eJeHHnqouavWr8mT\nJzcfMGBAm1tuueWywMDAuPnz54eWlJTgiSeeiIyOju7UpEmT2EGDBrU5ceKEJ3CuVeDFF18MDQ8P\n7xIWFtZl+vTpEc7ydmxBOHHihOfw4cNbh4eHdwkKCoq94YYbYgDg1KlTnn379m0bHBzcNSgoKLZv\n375t7U//fvDBB6N27NgROHXq1Jb+/v5xo0ePbgkAiYmJvr169bq8cePGsa1bt+707rvvBtvXm5qa\n6tmvX7+2gYGBcZ07d+6QlJRUYcvJ66+/HtK8efPOTZo0iX388ccjzfM2btzoHxsb295ms8WGhYV1\nGT16dMv8/HwBgO7du7cDgB49enT09/ePW7BgQXBF2+LMk08+GRkeHt4lICAgrnXr1p1Wr15tA4CK\nPgNn+8TPP//s06NHj3Y2my02ODi46+DBg9tUtM3OMJghIqqAn58fRowYgfj4+LJpK1asQPv27dG1\na1cAQEBAAOLj45Geno41a9bgzTffxKpVq8rls3nzZuzduxfr168/bx3h4eFISEhAZmYmFi1ahEmT\nJmHnzp1l81NTU5GRkYHjx49j4cKFuP/++3H2rH5g86OPPoodO3Zgy5YtSEtLwwsvvAAPDw8cP34c\ngwcPxrRp05Cenr5r7ty5x+68886Y5ORkr8zMTI8nnnii5Zo1a/bn5OQkbtmy5bcePXrkOtv+qVOn\nNtuxY0dAYmLir7t27fo1MTExYMqUKc26dOlSsGvXrl8AICMjI3Hr1q37HZfdtGnTPvv83NzcxBtu\nuCEHAHbv3h3Qrl27/LS0tF0PPfRQ6gMPPNC6tFQ/xPq2225r7eXlhaSkpJ8TExN/3bhxY+OXX37Z\nZTfNhg0bmgwfPvxsRkZG4r333nvmueeeC1+zZk2TTZs27UtJSdndpEmTkvHjx7d0+CxsBw4c+HnN\nmjW/z58/P9KdMSYjR468LC8vz+OXX3755dSpU7snTZp0AtAn7jFjxpw+cuTIT4cPH97j6+tbOmHC\nhJYAMH/+/OPdunXLnjNnzpHc3NzE+Pj4I5mZmR4DBw7808iRI9NOnz6964MPPkh67LHHWu7YscMX\nAMaPH9/K19e3NDk5efd777138MMPP3S57Tt27PB97LHHWi1cuPBgSkrK7jNnznidOHGikX2+l5cX\nXnrppaNpaWm7vv3229++/fZb2wsvvBAGANu3b98HAD/++OOvubm5iX//+9/PVrQtjnbv3u2zcOHC\n8G3btu3NyclJXL9+/f62bdsWAkBFn4GzfWLq1KnN+/Xrl5Genr7r+PHjex566KGTlX0ejhjMEBFV\nYsyYMVi5ciXy8/MBAPHx8RgzZkzZ/D59+qBz587w8PBAly5dcPvtt2Pz5s3l8pg5cyYCAgLg5+d3\nXv6DBw9GTEwMRAS9e/dG//79y7UEeXt7Y/r06fD29sagQYMQGBiIffv2obS0FO+99x5eeeUVREVF\nwdPTE7169YKPjw+WLl2KQYMGYdCgQfD09MTQoUMzO3XqlPPJJ580BgARUYmJiX7Z2dnSqlWrou7d\nu+c72/ZPPvkk5IknnkiJiooqbt68efG0adOSV65c2fRC6rN58+aFjzzyyGkvLy/cd999Z06dOuV9\n7Ngxr6NHj3pt3ry58TvvvHMkKCioNCoqqviBBx44sXLlSqetRgAQGxubc9ddd6V7enoiMDBQLVq0\nKGzWrFnHY2Jiivz8/NScOXOS165dG2weozFr1qyUoKCg0p49e+aNHDnyzLJly1zmDwCHDx/2/t//\n/td48eLFh8PCwkp8fHzU4MGDswEgMjKyZOzYsek2m600ODi49KmnnkrZtm2by+Bo+fLljaOiogoe\nfvjhM97e3rjmmmvyBg4cmL5s2bLg4uJirFu3rsns2bOTg4KCSnv06JE/YsSIM67y+vDDD4P79euX\nMXDgwGw/Pz81b968ZBFR9vl/+ctfcq+//vocb29vtGvXrnDs2LGnvvnmG5dlq8q2eHp6orCwUHbt\n2uVbUFAg7dq1K7ziiisKAMCdz8DMy8tLHTlyxOfQoUPe/v7+6sYbb8x2VUZXGMwQEVXi2muvRWho\nKFatWoWkpCRs27YNo0aNKpv/ww8/oG/fvggLC0Pjxo3x1ltv4fTp0+XyiI6Odpn/2rVrcdVVVyEk\nJARNmjTB559/Xm75pk2blhtn4+/vj+zsbJw+fRr5+fmIiYk5L8/Dhw/j448/RpMmTWCz2WJtNlvs\njh07AlNSUryDgoJKFy9e/Mc777wT1qxZs659+vRpm5iY6OusbKdOnWoUExNTYH/fpk2bwpMnT17Q\nwJSwsLCys5rNZisFgMzMTM8DBw40Ki4ulmbNmnW1l/mRRx5pdebMGZfra968eaH5fUpKSqM77rij\nrX35Tp06XeHp6Yljx46V5RETE1O2TKtWrQpSU1MboQJ//PGHd+PGjYvDwsLOu/wsKyvLY9SoUa2a\nN2/eOTAwMK5///7ts7KyPF11jR0+fLjRnj17Auzls9lssatWrQpJTU31Tk5O9iopKRHH8rkqV3Jy\nsndUVFRZ2qCgoNImTZqUrXjPnj0+ffv2bRsaGto1MDAwbs6cOVFpaWkuB2xVZVs6depU8Nxzzx19\n5plnmoeFhXW96aab2hw6dMgbcO8zMHvllVeOKaVw9dVXd2jbtu0V//nPf6ocLDOYISJyw+jRoxEf\nH4+lS5fixhtvRETEuaEWo0aNws0334yjR48iIyMD//jHP6CUKre8efCwWUFBAYYNG4ZHH30UJ06c\nQHp6OgYNGnTe8s6EhobC19cXSUlJ582Ljo7GXXfdhfT0dGRlZe3KysralZeXl/jcc8+lAsCwYcMy\nt2zZ8ntKSsruyy+/PH/8+PGtnK0jLCys0Dxu4+DBg43Cw8PduhTF1Ta70qZNm6JGjRqptLS0XfYy\nZ2dnJx44cOAXd9cRERFR9Omnn+63L5+VlbWroKBg52WXXVZW5qSkpLLg5ciRI40iIyPLBUTOypWR\nkeF1+vTp8y4/mzVrVsSBAwd8t27dujc7Ozvxiy+++A1A2ednbikBgOjo6KIePXpkmcuXm5ub+MEH\nHxxp3rx5saenp3Ion8sxM82aNSs6fvx4WdqsrCyP9PT0smBlwoQJrS6//PL833///afs7OzEqVOn\nHq9oOyvbFkf/+Mc/0nbs2LHv0KFDe0RETZw4sQVQ8WfgbJ9o2bJl8UcffXT45MmTe15//fXDjz/+\neKuqXmXFYIaIyA2jR4/Ghg0bsGDBgnJdTACQlZWFkJAQ+Pr6Ytu2bVi2bJnb+RYWFqKgoABhYWHw\n8vLC2rVr8cUXX7i1rIeHB8aNG4fJkycjOTkZJSUl+P7771FQUIA777wTn332GdavX4/i4mLk5uZK\nQkKCLSkpyfvo0aNeS5cubZKZmenh5+enAgMDSz08nJ8Ohg4dmjZ37txmycnJXikpKV6zZ89uNmzY\nMJddH2bNmjUr9vDwwN69e906MbVq1arommuuybj33nuj09LSPEpKSvDLL7/4rFmzJtCtCgFw9913\nn5w2bVqL/fv3NwKA5ORkr6VLlzYxp5kxY0azrKwsj+3bt/suX7489LbbbjtbWbmuu+66jLvvvrvl\nqVOnPAsKCmTt2rWBAJCVleXp6+tbGhoaWnLixAnPGTNmNDcvGxYWVvzHH3+Ubf+IESPSDx065Pv6\n66+HFBQUSEFBgWzevNl/586dvl5eXrjxxhvTn3rqqeZZWVkeO3bs8F2xYoXLVorbb7/97Ndff914\n/fr1gfn5+fLII480V0qVRQvZ2dmeQUFBJY0bNy5NTEz0fe+998LNyzdt2rR4//79ZWWrbFvMdu/e\n7fPf//7XlpeXJ/7+/srX11d5eHioyj4DZ/vEe++9F2wfaNy0adNiEYE9L3cxmCEickPr1q3Rq1cv\n5OTk4Oabby4374033sD06dNhs9kwa9YsjBgxwu18bTYbXn31VYwYMQLBwcFYtmzZeflX5MUXX0Tn\nzp3Ro0cPhISE4PHHH0dpaSmio6OxevVqPPfcc2jatGlsVFRUlxdffDGitLRUSktL5ZVXXomIiorq\n0qRJk9jvvvvO9tZbbx12lv/cuXNTunbtmtO1a9eOXbp06di5c+fcuXPnpri5baUPPvhgSu/evdvb\nbLbYr776KqCyZVasWHGosLBQOnTo0KlJkyaxw4cPjzl+/Ljb3VrTpk07OWjQoPT+/fv/KSAgIO7K\nK69sv3Xr1nLrve6667JiYmI6DRgwoN19992X+re//S2zsnyXL19+0NvbW7Vv375TWFhY15dffjkC\nAKZMmXIiPz/fIzQ0NPbKK6/s0L9//wzzchMnTjyRkJAQHBQUFDt27Njo4ODg0rVr1+7/+OOPQyIj\nI7tERER0ffzxx1vYrzJasGDBkZycHI9mzZp1HTt27GW33XbbaWflAYDu3bvnP//880fGjh17WWRk\nZNfg4ODiiIiIslamF1544egnn3wSEhgYGDd+/PhWQ4YMSTMv/69//St5woQJrW02W+y7774bXNm2\nmOXn53s8+eSTLUJDQ2MjIiK6nj592mvevHnHK/sMnO0T27ZtC7j66qs7+Pv7xw0dOrTts88+e6Rj\nx44VtpY5EneaMi8RrIgGorbvL+LKY3MmAAD+PfXtSlLWH3fvjWI1e/fuRYcOHeq7GBezHfVdgIZg\n3759jdq3b9+5sLBwB+9HYz27d+8O7dq1a2tn89gyQ0RERJbGYIaIiIgs7eK4pzYREVEl2rVrV6iU\nYpfbRYgtM0RERGRpDGaIqEHgxQhE5Irx/VDqaj6DGSKqd97e3mUPRyQiclRYWOgtIi4vFWcwQ0T1\nLjw8HMePH0dubi5baIionNLSUklOTg4sKSlZ7CoNBwATUb0LCgoCACQnJ8PVw+io+goLC10+eZnI\nAkoAfFtaWvqmqwQMZoioQQgKCioLaqjGta7vAhDVJnYzERERkaUxmCEiIiJLYzBDRERElsZghoiI\niCyNwQwRERFZGoMZIiIisjQGM0RERGRpvM8MkYXcs2R7tZZbOKZ7DZeEiKjhYMsMERERWRqDGSIi\nIrI0BjNERERkaQxmiIiIyNIYzBAREZGlMZghIiIiS2MwQ0RERJbGYIaIiIgsjcEMERERWRqDGSIi\nIrI0BjNERERkaQxmiIiIyNIYzBAREZGlMZghIiIiS2MwQ0RERJbGYIaIiIgsjcEMERERWRqDGSIi\nIrI0BjNERERkaQxmiIiIyNIYzBAREZGlMZghIiIiS2MwQ0RERJbGYIaIiIgsjcEMERERWRqDGSIi\nIrI0BjNERERkaQxmiIiIyNIYzBAREZGlMZghIiIiS2MwQ0RERJbGYIaIiIgsjcEMERERWRqDGSIi\nIrI0BjNERERkaQxmiIiIyNIYzBAREZGlMZghIiIiS2MwQ0RERJbGYIaIiIgszau+C0AXr3uWbK/v\nIhAR0SWALTNERERkaQxmiIiIyNIYzBAREZGlMZghIiIiS2MwQ0RERJbGYIaIiIgsjcEMERERWRqD\nGSIiIrI0BjNERERkaQxmiIiIyNIYzBAREZGlMZghIiIiS2MwQ0RERJbGYIaIiIgsjcEMERERWRqD\nGSIiIrI0BjNERERkaQxmiIiIyNIYzBAREZGlMZghIiIiS2MwQ0RERJbmVd8FIKLad8+S7dVabuGY\n7jVcEiKimseWGSIiIrI0BjNERERkaQxmiIiIyNIYzBAREZGlMZghIiIiS+PVTFSp6l4JQ0REVBfY\nMkNERESWxmCGiIiILI3BDBEREVkagxkiIiKyNAYzREREZGkMZoiIiMjSGMwQERGRpTGYISIiIkvj\nTfOIyKXq3jBx4ZjuNVwSIiLX2DJDRERElsZghoiIiCyNwQwRERFZGoMZIiIisjQGM0RERGRpDGaI\niIjI0hjMEBERkaUxmCEiIiJLYzBDRERElsZghoiIiCyNjzO4hFT31vREREQNGVtmiIiIyNIYzBAR\nEZGlMZghIiIiS2MwQ0RERJbGYIaIiIgsjcEMERERWRqDGSIiIrI0BjNERERkabxpHhHVuAu5QePC\nMd1rsCREdClgywwRERFZGoMZIiIisjQGM0RERGRpDGaIiIjI0hjMEBERkaUxmCEiIiJLYzBDRERE\nlsZghoiIiCyNN80jogalujfc4832iC5dbJkhIiIiS2MwQ0RERJbGYIaIiIgsjcEMERERWRqDGSIi\nIrI0BjNERERkaQxmiIiIyNIYzBAREZGlMZghIiIiS+MdgInoosA7BxNdutgyQ0RERJbGYIaIiIgs\njd1MRHRJY/cUkfWxZYaIiIgsjcEMERERWRqDGSIiIrI0BjNERERkaRwAbEHVHbBIRER0MWLLDBER\nEVkagxkiIiKyNHYzERHVId7XhqjmMZi5QBy/QkREVL/YzURERESWxmCGiIiILI3dTERE1cAuZqKG\ngy0zREREZGmilKrvMjQIIrIOQGgVFwsFcLoWilPTWM6axXLWLKuUE7BOWR3LeVopNaC+CkNU2xjM\nXAAR2a6UavDXS7KcNYvlrFlWKSdgnbJapZxENYXdTERERGRpDGaIiIjI0hjMXJh36rsAbmI5axbL\nWbOsUk7AOmW1SjmJagTHzBAREZGlsWWGiIiILI3BDBEREVkagxk3ichtIrJXRHJEJElE/mJMv15E\nfhORXBHZKCKt6rGMrUXkcxE5KyKpIvKaiHgZ82JFZIdRzh0iEluH5XpARLaLSIGILHaY57L+RMRH\nRN4TkUxjeybXRzlF5CoR+VJE0kTklIh8LCLNTPNFRJ4XkTPG63kRkboup0Oa6SKiROQG07QGUZ/G\nPH8ReUNETotIhoj8zzSvwdSniIwwjvssEflVRIY4zJ9k1GWmUbc+tVhOHxFZKCKHjfLsEpGBpvkN\n5lgiqmsMZtwgIn8F8DyAuwHYAFwH4A8RCQXwKYCnAIQA2A5geX2VE8AbAE4CaAYgFkBvAPeJSCMA\nqwEsBRAMYAmA1cb0upAM4FkA75knulF/MwFcDqAVgL4A/iUitXnjL6flhK6zdwC0NsqSBWCRaf69\nAIYA6AqgC4D/AzChHsoJABCRGAC3AkhxmDUTDaM+AV2fIQA6GH8nmeY1iPoUkSjoY2YygCAAjwFY\nJiLhxvwbAUwBcD10nbYB8HQtltMLwFHo47oxgGkAVhg/YhrasURUt5RSfFXyArAFwD1Opt8LYIvp\nfQCAPADt66mcewEMMr3/N4C3AfQHcBzGgG9j3hEAA+q4fM8CWOxu/UGfZPqb5j8D4KO6LqeT+X8G\nkOWwf9xren8PgK31VU4A6wAMAnAIwA2m6Q2iPgG0B5AJIMhF+gZRnwCuBHDSIc0pAFcb/y8D8Jxp\n3u8j/LkAAARMSURBVPUAUmu7nA7l2QNgWEM9lvjiq65ebJmphIh4AugOIExEDojIMaP7xg/AFQB2\n29MqpXIAJBnT68N/ANxmNOFHARgIfWK7AsAepZT50rU9qL9y2rmsPxEJhm5h2m1Kvxv1X2ZAt8z9\nYnpfbjtQj+UUkVsBFCilPneY3pDqsyeAwwCeNrqZfhKRYab5DaU+twPYKyI3i4in0cVUAH3sAM7L\nGSEiTeuicCISAeBP0PuiVY8lohrBYKZyEQC8AQwH8Bfo7ps46CbeQAAZDukzoLui6sP/oL+gMgEc\ng/4yXoWGV067isoVaHrvOK/eiEgXANOhuxzsHLcjA0BgbY7zcFE2G4DnADzsZHZDqs8WADoZ628O\n4AEAS0SkgzG/QdSnUqoEQDx0C0yB8XeCESi4KidQB3UqIt4APgCwRCn1m5Oy2MvTYI8loprEYKZy\necbf+UqpFKXUaQDzoJvxs6H70s2CoMdU1CkR8YBuhfkUuok5FHqsx/NoQOV0UFG5sk3vHefVCxFp\nC2AtgIeVUt+YZjluRxCAbIeWsLowE8D7SqlDTuY1pPrMA1AE4FmlVKFSajOAjdDdoUADqU9j8PQL\nAPoAaAQ9VuVdOTd43lk5gVquU+NYfx9AIXQg6Kws9vI0yGOJqKYxmKmEUuosdCuH+YvU/v8v0IMU\nAQAiEgAgBuW7IOpKCICWAF5TShUopc5AD1IdZJSni8Mv2y6on3Kauaw/o95TzPON/+ulzMaVIRsA\nPKOUet9hdrntQP2V83oADxlXq6QCiIYeIPp4A6vPPU6mmY+vhlKfsQD+p5TarpQqVUr9COAHAPYr\nxJyV84Rx7NUK4xheCN1iPEwpVeSsLA35WCKqDQxm3LMIwIMiEm70P08CkADg/wHoJCLDRMQXuvth\nj9HsW6eMFqODAP4pIl4i0gTAGOgTxyYAJdAnOh8Rsf+a+7ouymaUxxeAJwBPEfEVfcl4ZfUXD2Ca\niASLSHsAfwewuK7LaYw/+ho6UHzLyaLxACaLSJSINAfwSH2UEzqY6QR9Eo6FHvQ5AcDrpnLWe31C\nd4ceATDVSHMN9BU2603lbAj1+SOAv9hbYkQkDrqr2R6MxQO4R0Q6GsfbtNosp+FN6CvA/k8plWea\n3qCOJaI6V98jkK3wgh4z8waAdACpAF4F4GvMuwHAb9BN55sAtK7HcsYaZTgL4DSAFQAijHlxAHYY\n5dwJIK4OyzUT+pe3+TWzsvoD4AN9uWwmgBMAJtdHOQHMMP7PNr9Mywl0d0Sa8XoBpivH6rI+HdId\nQvmrmRpEfRrzrgDwPYAcAL8CGNoQ6xO6G+cAdJfMHwAecVh2slGXmdA/enxqsZytjLLlO+yLdxjz\nG8yxxBdfdf3is5mIiIjI0tjNRERERJbGYIaIiIgsjcEMERERWRqDGSIiIrI0BjNERERkaQxmiIiI\nyNIYzBAREZGlMZghIiIiS2MwQ0RERJb2/wHLk4HtApK53AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe04c162668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot\n",
    "plt.hist(\n",
    "    pp,\n",
    "    20,\n",
    "    label='Variances of the replicated data sets',\n",
    "    color=plot_tools.lighten('C0', 0.3)\n",
    ")\n",
    "plt.axvline(s2, color='red', label='Variance of the original data')\n",
    "plt.title(\n",
    "    'Light speed example with poorly chosen test statistic\\n'\n",
    "    r'$\\operatorname{Pr}(T(y_\\mathrm{rep},\\theta)\\leq T(y,\\theta)|y)=0.42$',\n",
    "    fontsize=16\n",
    ")\n",
    "plt.legend(loc='center left', bbox_to_anchor=(0.7, 0.8))\n",
    "plot_tools.modify_axes.only_x(plt.gca())\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
